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Abstract 

Asymptotic behavior of a class of nonlinear Schrodinger equations are studied. Par- 
ticular cases of ID weakly focusing and Bose-Einstein condensates are considered. 
A statistical approach is presented following [1] to describe the stationary probabil- 
ity density of a discretized finite system. Using a maximum entropy argument, the 
theory predicts that the statistical equilibrium is described by energy equivalued 
fluctuation modes around the coherent structure minimizing the Hamiltonian of 
the system. Good quantitative agreement is found with numerical simulations. In 
particular, the particle number spectral density follows an effective 1/k 2 law for the 
asymptotic large time averaged solutions. Transient dynamics from a given initial 
condition to the statistically steady regime shows rapid oscillation of the condensate. 

Resume 

Le comportement asymptotique des solutions d'equations differentielles hamiltoni- 
ennes est presente dans le cas general des equations de Schrodinger nonlineaires. Ce 
travail reprend une etude precedente s'appuyant sur une description statistique de 
l'espace des phases de la solution [1]. La recherche de la distribution stationnaire a 
l'equilibre statistique s'effectue pour la dynamique discrete en maximisant l'entropie 
autour de la solution concentrant toute la masse du systeme. On trouve alors que 
la distribution d'equilibre correspond a l'equipartition statistique de l'energie en 
exces sur tous les modes accessibles. Les simulations numeriques sur un modeles 
faiblement focalisant et dans le cas particulier d'un modele ID de condensat de 
Bose-Einstein permettent de montrer un bon accord quantitatif avec les predictions 
de la theorie. 

Key words: Wave turbulence, Bose-Einstein condensate, statistical equilibrium 
Mots-cles : turbulence d'ondes, condensats, distribution d'equilibre 
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Version francaise abregee 
1 Introduction 

The emergence and persistence of large scale coherent structures in the midst 
of small scale turbulent fluctuations is a common feature of many turbu- 
lent fluids, plasma systems and nonlinear dynamics in general [2,3,4,5]. Clas- 
sical examples include the formation of large vortices in two dimensional 
high Reynolds turbulence or the emergence of solitary solutions in nonlinear 
optics[6,7]. Moreover and quite surprisingly, such complex phenomena survive 
even when the dissipationless versions of the dynamics are taken. In these 
cases as in many classical Hamiltonian systems, the dynamics are formally 
reversible [8]. Therefore one can question how such apparently irreversible pro- 
cesses can be compatible with these irreversible dynamics? The recent at- 
tainment of Bose-Einstein condensates (BEC) [9,10,11] provides a fascinating 
new system exhibiting similar behavior. After being attained, the condensate 
can indeed be viewed as a large coherent structure, persisting in the atomic 
trap and surrounded by fluctuations. When thermal fluctuations are neglected, 
BEC can be modelled by the usual semi-classical Gross-Pitaevskii equation. 
It is strictly valid only at T = and is also Hamiltonian dynamics. Moreover 
note that for all these examples the dynamics resume in partial differential 
equations describing the evolution of one or few classical fields (velocity for 
turbulent flows or wave function for BEC for instance). 

The interplay between the fluctuations and the coherent structures in these 
systems is of crucial interest to our understanding of nonlinear dynamics. Im- 
portant questions are related to the asymptotic behavior of the dynamics and 
to their possible statistical description. In particular, the expected thermal 
equilibrium would in fact lead to the well-known Rayleigh- Jeans divergences 
for a classical field. The aim of this note is precisely to investigate how this 
black body-like catastrophy manifests in the a priori smooth dynamics driven 
by the partial differential equations (PDE) considered here. For simplified re- 
versible models, it is believed and borne out by numerical simulations that 
the coherent structures may act as statistical attractors to which the whole 
system relaxes. Following weak turbulence theory, in continuum systems, the 
fluctuations were shown to decrease without bound in addition to cascading 
to smaller and smaller scales. The explanations proposed for this cascade-like 
process invoke in fact thermodynamical considerations[12,13]. In collabora- 
tion with R. Jordan, we have recently focused on such scenario for the class 
of nonintegrable nonlinear Schrodinger (NLS) equationsfl]. Following [21] we 



Email address: josseran@lmm.jussieu.fr (Christophe Josserand). 



2 



seek the probability density of a solution for a finite number n of modes version 
of the NLS equations. Usual thermodynamical arguments, where the probabil- 
ity density corresponds to the standard Maxwell-Boltzmann distribution, fail 
and one has to consider that a coherent structure emerges from the solution 
of the PDE. Indeed, from numerical simulations we assume the presence of 
a non-zero mean field which contains most of the conserved particle number 
(L 2 -norm squared). Such specific treatment of a particular part of the solution 
is therefore very close to the statistical theory of the Bose-Einstein transition! 
The stationary probability distribution is obtained via the maximization of 
the entropy of this finite statistical system. It describes an ensemble where 
the mean-field corresponds to a large-scale coherent solitary wave, which min- 
imizes the Hamiltonian given the particle number, coupled with small-scale 
random fluctuations, or radiation. The fluctuations equally share the differ- 
ence of the conserved value of the Hamiltonian and the Hamiltonian of the 
coherent state. The effective temperature of this thermal-like system is in- 
versely proportionnal to n the number of modes and goes therefore to zero 
in the continuum limit. Thus the discretization level n of the dynamics trig- 
gers an effective energy cut-off which avoids the Rayleigh- Jeans divergencies. 
This statistical theory found very good qualitative and quantitative agreement 
with numerical simulations done for a weakly focusing NLS equation. Statis- 
tical ensemble can be retrieved by time and ensemble averaging of large time 
numerical solutions starting from random initial conditions. At large enough 
time, the dynamics of the solutions present statistical stationary properties 
in full agreement with the statistical theory. The road to this statistically 
stationary state is also investigated and allows a consistent scenario for the 
dynamics in the continuum limit n — > oo. 

In this note, I discuss how these results apply to the case of BEC dynamics. 
The Gross-Pitaevski equation in harmonic trap used to model the condensate 
evolution is in fact a particular case of NLS system so that the introduced 
description is valid. Before that, I will present in detail the general theory for 
the NLS equation in one spatial dimension with a finite number of modes. 



2 Self-organization in NLS-systems 



In this section I briefly introduce the results obtained for the NLS-equations. 
Most of the conclusions and the figures are presented for a specific 1-D focusing 
NLS equation but the results apply to the whole class of these models. The 
presentation here will follow [1] and the reader will find more details and 
references there. 
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2. 1 Generalities 



We consider the general dimensionless NLS equation: 



id^ + A^ + f(\^ = 0, 



(1) 



where ^(r, t) is a complex field and A is the Laplacian operator. The function 
f{\ip\ 2 ) stands for nonlinear interaction and external potential. Among other 
phenomena, it is used to model gravity waves on deep water [14], Langmuir 
waves in plasmas [7], pulse propagation along optical fibers [4], superfluid 
dynamics[15] and Bose-Einstein condensates [16]. In this latter case, the inter- 
action function f{\ip\ 2 ) depends also on the position r to describe the atomic 
trap in which Bose-Einstein condensates are achieved. When f(\ip\ 2 ) = ±|"0| 2 
and eqn. (1) is posed on the whole real line or on a bounded interval with 
periodic boundary conditions, the equation is completely integrable [18]. In 
any other configuration it is nonintegrable. 

The Hamiltonian equation associated to (1) is: idtip = 5H/5ip*, where ip* is 
the complex conjugate of the field ■0, and H is the Hamiltonian: 



Here, in addition to the kinetic term (V^l 2 , the potential F is defined via the 
relation F(a) = Jq f(y)dy. The dynamics (1) conserves, in addition to the 
Hamiltonian, the particle number (also called the mass) 



Without loss of generality we hereafter restrict the statistical analysis and the 
numerics to non-integrable models in one spatial dimension. Equations (1) 
exhibit solitary wave solutions ip = (fi(x,t)e a 1 which satisfy 



It is worthwhile to notice that the solution of equation (4) minimizes the 
Hamiltonian (2) for a given particle number. We will denote H*(N) the value 
of the energy of this solution. These localized structures are found to play 
an important role in the time evolution of Eq. (1). Figure (1) shows indeed 
snapshot of the dynamics for the weakly focusing nonlinearity f{\ip\ 2 ) = 
Starting with a slightly perturbed homogenous solution, a collection of solitary 
peaks emerge from this linearly unstable state. These solitary waves rapidly 




(2) 




(3) 



0^ + /(|0| 2 )0-A 2 = O. 



(4) 
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Figure 1. Profile of the modulus |^| 2 at four different times for the system (1) 
with nonlinearity f(\ip\ 2 ) = \ip\ and periodic boundary conditions on the interval 
[0,256]. The initial condition is ip(x,t = 0) = A, with A = 0.5, plus a small ran- 
dom perturbation. The numerical scheme used to approximate the solution is the 
split-step Fourier method. The grid size is dx = 0.125, and the number of modes 
is n = 2048. a) t = 50 unit time: Due to the modulational instability, an array 
of soliton-like structures separated by the typical distance ti = 2tt/\/ A/2 = 4ir is 
created; b) t = 1050 unit time: The solitons interact and coalesce, giving rise to a 
smaller number of solitons of larger amplitude; c) t = 15050: The coarsening process 
has ended. One large soliton remains in a background of small-amplitude radiation. 
Notice that for t = 55050 unit time (Figure d)), the amplitude of the fluctuations 
has diminished while the amplitude of the soliton has increased. 

coalesce into a single coherent structure surrounded by small amplitude fluc- 
tuations. A slow, quasistatic dynamics follows where the fluctuations are seen 
to decrease in wavelength and in amplitude. The solitary structure gathers 
almost all the mass of the system while the energy reaches smaller scales. 
In other words, the dynamics presents the condensation of the mass into a 
coherent structure while the energy is distributed in the system. 



2.2 Statistical model 



To explain these numerical observations, we construct a consistent mean-field 
statistical theory based on [21]. It seeks to describe the solution as a mean field 
which contains most of the particle number of the system while the amount 
of energy not accounted by this solution is dispersed towards small scales. 
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The theory is in fact presented for finite dimensional approximation of the 
NLS equation, using also Dirichlet boundary conditions on the interval [0,L]. 
Other boundary conditions, such as the periodic ones used in the simulations 
can easily be considered without any changes in the conclusions. 

Let ej(x) = \JljL €m{kjX) with kj = nj/L, and for any function g(x) on f2 = 
[0, L] denote by gj = J Q g(x)ej(x) dx its jth Fourier coefficient with respect to 
the orthonormal basis ej,j — 1, 2, • • •. We consider now a truncated version of 
the NLS equations to the n first Dirichlet modes. We write tp^ = + iv^ n \ 
where and v^> are two real functions whose Fourier coefficients satisfy 
the coupled system of ordinary differential equations 

% - fyj + (f((u (n) ) 2 + (u (n) )> (n) ) = (5) 
vj + k)u 3 - (f((u^) 2 + (v^) 2 ^)^ = . (6) 

It corresponds for tp^ to: 

<W B) + ^£ ) + ^(W (B) lV B) ) = o, 



where P n is the projection onto the span of the eigenfunctions e±, ■ ■ • , e n . This 
equation is a natural spectral approximation of the NLS equation (1), and it 
may be shown that its solutions converge as n — > oo to solutions of (1) [17,19]. 

For given n, the system of equations (6) defines a dynamics on the 2n- 
dimensional phase space R 2ri . This finite-dimensional dynamical system is a 
Hamiltonian system, with conjugate variables Uj and Vj, and with Hamiltonian 

H n = K n + Q n , (7) 
where 

1 r 1 n 

K n = -j ((4 n) ) 2 + (^ n) ) 2 ) dx = - £ k 2 (u 2 + v 2 ) , (8) 

is the kinetic energy, and 

Q n = -± J F((uW) 2 + (vW) 2 )dx , (9) 

is the potential energy. The Hamiltonian H n is, of course, an invariant of the 
dynamics. The truncated version of the particle number, up to a multiplicative 
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factor in the definition 



1 r 1 n 

n = ~ J ((U^r + (V^f) dx = t £( U J + „J) , (10) 



is also conserved by the dynamics (6). 

We argue that we can build a statistical treatment of the system using the 
usual assumption that the dynamics is ergodic and noting also that system 
(6) satisfies the Liouville property (the measure Y¥- =1 dujdvj is invariant) [20]. 
We thus replace the time dependant dynamics by a statistical (time inde- 
pendant) description of the solution at any time. We introduce a probability 
density p^ n \ui, • • • , u n , V\ • • • , v n ) on the 2n-dimensional phase space and we 
seek the density function p^ which maximizes the Gibbs-Boltzmann entropy 
functional: 

„ n 

S(p) = - J plogpYldujdvj , (11) 

R 2n 3 = 1 



We easily obtain the usual canonical ensemble solution: 
p oc exp (-(3H n - pN n ) , 



subject to the mean constraints (H n ) = H° and (N n ) = N°. H° and 
are the given values of the Hamiltonian and the particle number, respec- 
tively, and (3 and p are the Lagrange multipliers associated to these con- 
straints. Such a density is actually ill-defined since it is not normalizable (i.e., 
J R 2n exp[— (3H n — pN n ] Ilj=i dujdvj diverges) for the focusing nonlinearities we 
discuss here [21, 22]. Moreover, this distribution does not take into account the 
numerical observed fact that an important part of the phase space consists 
of configurations where most of the mass is concentrated in the solitary wave 
solution. 

Inspired by these remarks, we have built an adapted statistical description 
of the dynamics. We decompose the fields u n and v n into two contributions: 
the means denoted (uj) and (vj) and the fluctuations (5u^ n \5v^) = (u^ — 
-i/™) — (*/"))). This decomposition will be clarified by the yet to be de- 
termined ensemble p^ n \ We now impose that in the long-time NLS dynamics 
(6) and for the continuum limit n — > oo, the amplitude of the fluctuations van- 
ishes. Consequently, the number of particles in this limit is almost determined 
by the mean field. The vanishing fluctuations hypothesis is written: 

n 

((5u^) 2 ) + ((5v^) 2 )} dx = J2 \((Suj) 2 ) + ((Svj) 2 )} -> , as n -> oo .(12) 
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Using this assumption, it follows that for n large enough, the total number 
of particles is well approximated by the mean contribution. Similarly, one can 
show that the potential energy is almost entirely determined by the potential 
of the mean. However, although the fluctuations do not contribute much to the 
particle number and the potential energy, they may have a significant kinetic 
energy. Indeed, this contribution (1/2) £"=i k 2 [{(Suj) 2 ) + ((5vj) 2 )], need not 
tend to as n — > oo, even if (12) holds. Actually, it must in general not tend 
to since the total energy of the system is conserved. 

From these remarks, we can impose the following mean-field constraints on p: 

N n (p) = ^E((^) 2 + (v 3 ) 2 ) = N° (13) 

1 n 1 r 

Hn{ P ) = ^ E k Wj) + <«;» - 2 / + ^ (n) > 2 ) dx = H ° ■ ( 14 ) 

Where iV and H° are precisely the conserved values of the particle number 
and the energy, as determined from initial conditions. 

The solutions of the statistical equilibrium states with the particle number 
constraint (13) are calculated now invoking again a maximum entropy princi- 
ple. Notice that, similar to its above use, this principle has no reason to hold 
in such hamiltonian and reversible systems. However, it allows a consistent 
calculation of the density following an ergodic assumption [20]. It somehow 
corresponds to the determination of the density distribution around the most 
probable state. 

The solutions p*™-* are therefore calculated through usual techniques: two La- 
grange multipliers are used to enforce the mass and the energy constraints. 
Considering independent statistical variables, the factorization of the maxi- 
mum entropy distribution is straightforward: 

n 

p {n) (u 1 ,...,u n ,v 1 ,...,v n ) = JJpjiuj^j), (15) 

j'=i 

where, for j = 1, . . . , n, 

Pj(uj, vj) = ^ exp j-^ 1 ((uj - (uj)) 2 + (vj - ( Vj )) 2 ) | , (16) 

In addition, we find that the complex field {ip^) = (u^)+i(v^) is a solution 
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of (setting A = fi/[3) 



(17) 



which is clearly the spectral truncation of the eigenvalue equation (4) for the 
continuous NLS system (1). It follows, therefore, that the mean-field predicted 
by our theory corresponds to a solitary wave solution of the NLS equation. 
On the other hand, the fluctuations 8uj and 5vj are independant Gaussian 
variables with identical variance 4ri- 

The energy constraint (13) imposes: 

H° = ^ + H n ((u^),(v^)). (18) 



where the term ^ reflects the equipartition of energy among the 2n indepen- 
dant fluctuating modes. The calculation of the total entropy is straightforward: 

= C(n) + n log ( L2 ' g0 - g "<<"'"'-<"" ,, »l) . (19) 



where C{n) = n — YTj=\ l°g(j ' 7r /2) depends only on the number of Fourier 
modes n. Thus, the key results follow the maximization of the entropy: firstly 
the mean-field pair ((w^), (?/"))) has in fact to realize the minimum possible 
value of H n over all fields (u^ n \ v *")) that satisfy the constraint N n (u^ n \ v^) = 
N°. Moreover the excess energy not present in the mean field is equally dis- 
tributed among the fluctuating modes, with the inverse temperature: 

o-n^m- (20) 



where H* is precisely this mimimum value of H n allowed by the particle num- 
ber constraint N n = N°. We obtain also that the inverse temperature scales 
linearly with the number of modes in the continuum limit. 
The vanishing of fluctuations hypothesis is also verified if one computes the 
contribution of the fluctuations to the number of particles. It reads: 

i n ttO tt* "1 1 

^£[<(^) 2 > + <(^) 2 >] = a Ep = °(-). asn^oo. (21) 



The limitation of our statistical description follows directly from Eq. (21). In- 
deed, for a fixed number of modes n, the coherent structure can only emerge 
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if N 3> — — — . If this criterion is not satisfied, the condensation of the solu- 
tion into the coherent structure can even be broken. In fact, it is still under 
debate whether such a coherent structure can emerge in the continuum limit 
when starting with a highly fluctuating state concentrated at relatively small 
wavelengths [24]. In this approach the fluctuation modes have been choosen as 
the ej[x) Fourier functions. This assumption is in fact a direct consequence of 
the factorization of the entropy, considering independent Fourier coefficients. 
Such an hypothesis is correct for defocusing homogenous NLS equation but 
is in general only valid for large wave number kj. Indeed, it is well known by 
linearizing eq. (1) around a given state that the perturbation modes have to 
be found through the so called Bogoliubov theory. Although the full determi- 
nation of the Bogoliubov modes is needed to improve the statistics, we note 
that plane waves are a good approximation of the fluctuation modes for large 
enough kj so that our theory is always relevant to describe the statistics of 
the small wavelength perturbations. 

Finally, regardless of these restrictions, our statistical equilibrium approach 
gives the following prediction for the particle number spectral density: 

m 2 ) = m\ 2 +^0^i (22) 

where we have used the identity ipj = Uj + ivj. 



2.3 Numerical results 

In [1], we discuss how numerical simulations compare with these statistical 
predictions. Starting from any set of initial condition, one expects, follow- 
ing ergodic theory, that large time and large ensemble (over different initial 
conditions) averages will reproduce the statistical description of the system. 
Alternatively, one can choose to measure "quasi-static" averaged quantities 
by measuring only few time unit averages and (if possible) large ensemble 
average. Such quasi-static description would converge to the sought after sta- 
tistical regime for large enough integration time of the dynamics where the 
transitory effects due to the initial conditions can be ignored. We address here 
the weakly focusing nonlinearity f{\i>\ 2 ) = \ip\ already used for figure (1). The 
transient dynamics can be estimated by plotting the two contributions to the 
conserved total energy. Figure (2) shows over a large period the evolution of the 
kinetic and the potential energy, averaged over 16 slightly noisy homogenous 
states. Fast but small oscillations of the energies around a smooth evolution 
account for the rapid fluctuation modes. The kinetic (potential) energy is then 
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Figure 2. Kinetic (top curve) and potential (bottom curve) energies as function 
of time for the weakly focusing NLS. The straight line below the potential energy 
corresponds to the potential energy of the minimizer of the Hamiltonian for the 
particle number. The calculation are made with n = 512 discretization points and 
ensemble averaged is performed over 16 different initial configurations. 

found to increase (decrease) as time goes on, indicating the transfer of energy 
towards smaller scale. Moreover, after a short transient (until t = 50000 time 
units) corresponding to the formation of peaks and to their coalescence into 
one single structure, a slow quasistatic regime is observed as described in figure 
(1). This slow evolution witnesses the convergence of the coherent structure 
to the Hamiltonian minimizer. The straight line below the potential energy 
indicates the potential energy of the minimizer for the total particle number 
No. The potential energy approaches a constant value close to this bound, 
the difference being due to the finite particle number which still remains in 
the fluctuations. This quasistatic dynamics seems in fact to saturate for this 
n- modes calculations at large time (t > 900000). Indeed, such convergence is 
observed on the mean particle number spectra averaged over few time units 
and 16 initial conditions, which shows stationary properties thereafter. 

Figure (3) shows this particle number spectral density. The spectral density of 
the solitary solution which contains the total initial particle number is drawn 
(smooth line) for comparison. The equipartition for large wavenumber deduced 
from the theory is also marked, corresponding to a straight line in the log-log 
plot. Very good quantitative agreements are found between the theoretical 
predictions and the numerics, for the coherent structure at large scale as well 
as for the 1/k 2 spectrum for small scale. 

Before this stationary regime is attained, one can ask how the system finds 
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Wave number k 



Figure 3. Particle number spectral density \tpk\ 2 as a function of k for t = 1.1 x 10 6 
unit time (upper curve). The lower curve (smooth one) is the particle number spec- 
tral density for the solitary wave that contains all the particles of the system. The 
straight line drawn for large k corresponds to the statistical prediction (22) for the 
spectral density for large wavenumbers. The numerical simulation has been per- 
formed with n = 512, dx = 0.25, N° = 20.48 and H° = -5.46. 

its road to the statistical equilibrium starting from a given initial conditions 
(taken here as homogenous states perturbed by a slight random noise). Since 
the dynamics after the coalescence regime is quasistatic, converging slowly 
towards the statistical equilibrium, one can follow this evolution using the av- 
erage of the particle number spectral density at time t. Here the mean value 
is taken over the ensemble of 16 initial conditions and also by time averaging 
over few time units only. Thus the average smoothes the fast variations while 
the evolutions due to the quasistatic dynamics can be neglected. Figure (4) 
shows the particle spectrum for such an intermediate time tj. The statistical 
equilibrium has not been reached yet, but the solution has already converged 
into a single solitary wave surrounded by fluctuations. The spectrum at long 
wavelengths accounts again for the soliton-like coherent structure containing 
already almost all the total particles. On the other hand at small scales, the 
fluctuation modes already exhibit the 1/k 2 law but only for wavenumbers 
smaller than a well defined wavenumber K(ti). For higher wavenumbers, the 
amplitudes of the fluctuation modes are at the level of the initial noise. Some- 
how, the system acts as if at each time t the statistical equilibrium (22) were 
satisfied only for scales larger than X(t) = 2n/K(t). As time goes on, the front 
k = K(t) advances in the momentum space as the "quasistatic" equilibrium 
invades smaller and smaller lengthscales. For a finite number of modes, this 
slow dynamics saturates when all avalaible modes are reached by the equi- 
librium. Thermal equilibrium between the fluctuations has ben attained. The 
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Figure 4. The particle number spectral density for n = 512 and dx = 0.1 (thus 
L = 51.2) at unit time t = 5 ■ 10 5 . The coherent soliton structure already accounts 
for almost the entire number of particles of the system, but the system has not 
yet reached statistical equilibrium. The initial noise level is still present for large 
wavenumbers (k > 20), while at smaller wavenumbers, one can recognize both the 
soliton-like structure and a fluctuation spectrum following approximately a k~ 2 law. 
The spectrum has been obtained by an ensemble average over 16 initial conditions 
and a time average over the final 10 unit times. 

front dynamics is found to obey the scaling K(t) ~ t 1 / 4 using high derivative 
moments of the wavefunction ip. It is then tempting to generalize these results 
to the continuum limit n — > oo. In that case, we expect that at any time the 
dynamics is made of a solitary coherent structure containing most of the parti- 
cle number in quasi-statistical equilibrium with a finite number of fluctuating 
modes. As time goes on, the number of modes at equilibrium increases and 
the effective temperature due to the equipartition of the energy among the 
modes decreases. Thus the contribution of the fluctuations to the total mass 
decreases providing an inverse transfer of mass from the fluctuations to the 
soliton-like structure. The t — > oo limit corresponds to the coherent structure 
containing all the mass of the system while the excess energy is distributed 
over an infinite number of modes of zero amplitude! 

One may have argued that the prediction of the equipartition of the energy 
among the modes surrounding a condensate structures was obvious since some 
sort of thermal equilibrium was assumed. In fact, it is noteworthy to observe 
that the dynamics of (1) reaches in fact a self-thermalized state through the 
dependance of the temperature with the number n of available modes. The 
singular limit n — > oo can then be interpreted as a semi-classical example of 
the Rayleigh- Jeans paradox[13]. 
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3 Wave turbulence in BEC 



We outline here the conclusions of the previous section for a model of BEC. For 
sake of simplicity and to allow long numerical simulations we will again restrict 
our study to one spatial dimension. In dimensionless units, the dynamics of a 
BEC trapped in a harmonic potential reads[16]: 

icki>= (-^A + V ext (x) + M 2 )i> (23) 



where the additional harmonic potential V ext (x) = ^Q 2 x 2 describes the exter- 
nal potential used to confine the particles in an atomic trap. The usual de- 
focusing nonlinearity corresponding to repulsive particle interactionsis used. 
The integrability of the ID NLS equation is broken here due to the presence of 
the external potential V ext . The structure of the equation is clearly equivalent 
to the general NLS system. The number of particles: 



N = J \^\ 2 dx 
is conserved as well as the Hamiltonian 



H=±J (|VVf + fiV|Vf + H 4 )^. (24) 



To the kinetic energy |Vt/>| 2 and the nonlinear potential term \ip\ A is added 
the contribution of the external potential fl 2 x 2 \i(j\ 2 . The solitary solutions 
<j)(x)e~ lflt which minimize the Hamiltonian for a given number of particles N 
are solutions of: 

H> = -\<t>*x + \&x 2 4> + 3 = . (25) 



Such solutions can be obtained numerically and are well described for large 
enough N by the so-called Thomas-Fermi approximation. Neglecting the ki- 
netic term the Thomas- Fermi solution i^tf only balances the nonlinearity and 
the external potential: 



for |x| < and 4>tf{x) = elsewhere. By a straightforward integration the 
chemical potential /itf satisfies: 
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Figure 5. Solutions <j){x) of Eq. (25) for 17 = 0.1 and different particle numbers 
N = 11.5, 45, 180 and 720 (as expected, the smaller N the smaller <j> on the figure). 
For each solution the dashed line show the Thomas-Fermi approximation 4>tf- 



where N is the total particle number. The radius Rtf of the condensate in 
the Thomas-Fermi approximation is found: Rtf — {3N/(2Q 2 )) 1 ^ 3 . Figure (5) 
shows solitary solutions of Eq. (25) for different values of iV together with the 
corresponding Thomas-Fermi solution. We observe that the Thomas-Fermi 
approximation is very good except in a small boundary layer near x = Rtf- 
The potential \i is also found numerically very close to iitf- The solutions were 
calculated on the periodic box x £ [—64, 64] of length L = 128 unit length, 
with O = 0.1, that we will keep constant throughall later on. 

Perturbations around this ground state are studied to determine the eigen- 
frequencies of the fluctuation modes. Following Bogoliubov approach, we are 
seeking solutions for the weak perturbation regime: 



At first order in the complex functions u and v, we obtain from (23) the 
coupled equations: 




V>(x,t) = (4>{x) + u(x)e~ iujt + v(x)e iujt ) e 



—ifj,t 
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[~2 dxx + ~ ^ + 2 ^ 2 J U ^ + ^ x f v ^ x ) ( 26 ) 

— ouv(x) = {-\d xx + -Q 2 x 2 -/! + 2<p(x) 2 ) v(x) + <j>(x) 2 u(x) (27) 



This linear system can be solved numerically to determine the excitation 
spectrum[23]. However, as for the NLS system described above, high frequency 
modes are well approximated by plane waves u(x), v(x) oc e lkx for small wave- 
lengths such that k 2 ^> fi ~ Htf- We retrieve then the well-known Bogoliubov 
relation for high wavenumber dispersive waves 




A priori, the statistical approach developped above applies to this specific 
NLS system (23). Thus, starting from any initial condition, we should observe 
the formation of a coherent structure solution of (25) containing almost all 
the particle number of the system, in the midst of wave fluctuations. For large 
enough time, one should observe also the statistical equipartition of the excess 
energy between all the avalaible modes for a finite-n grid point calculation. 
In particular, we should observe again a 1/k 2 particle number spectrum for 
large wavenumbers. To illustrate this dynamics, we initiate the dynamics with 
the ground state of a condensate of noninteracting bosons for N = 180. When 
the nonlinear term is neglected, the ground state of the linear Schrodinger 
equation is described by a gaussian distribution of the particle around the 
trap center: 




The gaussian width a = 1/ defining the typical radius of the noninteracting 
gas. Numerical simulations are typically made on a regular grid of 512 or 1024 
points, using a Yoshida pseudospectral splitting scheme [25] which is 4th order 
in time. Mass and energy are conserved to within 10 -10 and 10~ 5 relative 
errors respectively. At t = we start the dynamics of Eq. (23) with: 



ip(x, 0) = ipi(x, 0) + r)(x) 

where r](x) is a white noise of very small amplitude (10~ 5 ), taken to break the 
x — > —a; symetry of the system. Figure (6) show snapshots of the solution for 
different times. 
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Figure 6. Particle number density \tp(x,t)\ 2 for different times: a) initial gaussian 
profile at t = 0, b) t = 4000, c) t = 34000 and d) t = 250000 units time. The 
coherent structures containing all the particles is shown as dashed lines on figure 
d). 



We observe roughly the expected dynamics where the coherent structure 
emerges from fluctuations. However, one can notice that the fluctuation am- 
plitudes are much higher here than for the simulations shown in the previous 
section, so that the comparison between the solution at large times and the 
coherent structure (figure 6 d)) is only partially satisfied. This high level of 
fluctuations is due to the large difference of energy between the initial condi- 
tion (H = 564.25) and the coherent structure <f>(x) minimizing the Hamilto- 
nian (H* = 486.15). Moreover, since the initial gaussian width (<r = 3.16) is 
smaller than the Thomas-Fermi radius (Rtf — 30) of the coherent structure, 
we observe in the numerics numerous large radius oscillations around Rtf 
before the solution stabilizes. Such radius oscilations can be noticed between 
figures (6 a), b), c)). 

These oscillations and their "effective" damping are even better seen on figure 
(7) where the different (kinetic, external potential and nonlinear) contribu- 
tions to the total energy are shown at short time. Recall that the sum of 
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Figure 7. Different contributions to the total constant energy as functions of time for 
short times. Kinetic energy (black curve), external potential contribution (red curve) 
and nonlinear one (blue curve) have been also vertically translated for vizualisation. 
The calculation is made with n = 512 discretization points. 

these three contributions is constant throughout the dynamics. The ampli- 
tude of these oscillations are rapidly decreasing after few thousand time units. 
Thereafter, only small oscillations are observed around slowly varying energy 
contributions. This slow dynamics is shown on figure (8) where again kinetic 
(bottom curve), and the sum of the external and the potential energies (top 
curve) are drawn on a much larger time scale. After the short transient, where 
fluctuations are important, we observe a quasistatic dynamics consisting of 
small rapid fluctuations around a slow variation of the energy contribution. 
As expected, the kinetic contribution increases, while obviously the other con- 
tributions decrease. This slow dynamics indicates as for the previous focusing 
NLS system the transfer of energy towards smaller and smaller scales as time 
goes on. 

The cascade-like transfer can also be seen on figure (9) where the particle 
number spectral density of the fluctuations are shown at different instants of 
the dynamics. This spectral density is obtained by substracting the coherent 
structure <pe~ lflt from the solution ijj(x,t). The spectra are obtained by time 
averaging over one unit of time, except for the initial condition spectrum. One 
can see the evolution of the spectrum from the initial gaussian distribution 
to the predicted spectrum at t — 422000. The dashed line represents the 
quantitative prediction (\ipj\ 2 ) = - ~ k ^ n , based on H* calculated with the 

i 

soliton-like solution containing all the particles (N = 180). Such a stationary 
spectrum corresponding to the statistical equilibrium distribution was actually 
mostly attained since t = 200000 time units. The curve at t = 16900 indicates 
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Figure 8. Kinetic energy (bottom black curve), sum of the external potential and 
of the nonlinear contribution (top red curve) as function of time. The calculations 
are made here with n = 1024 discretization points. 

that the 1/k 2 spectrum is fist reached at long wavelengths and invades smaller 
and smaller lengthscales as time goes on, as already observed for the focusing 
NLS of the previous section. 

A difference between the predicted spectrum amplitude at high wave numbers 
and the numerical results is observed, while the agreement for the focusing 
NLS in the previous section was much better. We understand such discrep- 
ancy by a higher level of fluctuations in the BEC system. The energy difference 
between the initial condition and the ground state coherent structure contain- 
ing all the particles is AH = 78 here while it was only a few units for the 
focusing NLS case. Thus the amplitude of the fluctuations are the same order 
as the coherent structure, as observed on figure (10) where, both the density 
of the numerical simulations and the density of the fluctuation field deduced 
after sustraction of the expected coherent structures are shown at t — 422000. 
Here we are close to the limitation of our approach since the particle num- 
ber of the fluctuations field is close to the total particle number as discussed 
previously. 

This finite mode self-thermalization of the NLS-like equations has been used 
recently to model finite temperature BEC[26, 27,28]. Distributing fluctuation 
energy equally among the Bogoliubov modes, the authors "mimic" numerically 
a thermal system at equilibrium with only the Gross-PitaevskiT equation. Al- 
though the cutoff at high wavenumbers of the black-body-like catastrophe is 
needed, we want to emphasize here that this approach is inconsistent since 
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Figure 9. Particle number spectral density of the fluctuations around the solution 
(p(x) e —ifit at different stage of the dynamics on a log-log scale. From bottom to top 
at high wave number, we have the spectrum of the initial condition t = (black 
curve), the spectrum at t = 16900 time units (red curve) and at t = 422000 (blue 
curve). The dashed line represents the spectral density for large wavenumber as 
predicted by our theory. The number of grid points for calculation is n = 1024 here. 
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Figure 10. Particle density \ijj(x,t)\ 2 (solid line) of the solution of the dynamics (23) 
and the density of the fluctuation field (dashed curve) at t = 422000. 
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no convergence of this regularization in the continuum limit is proposed. A 
consistent regularization would have to integrate a quantum description of the 
dynamics. Our work here, without solving this crucial question, can be seen 
as an attempt to describe the dynamics of a classical field. 



4 Conclusion 

We have presented here the asymptotic dynamics of different nonlinear dif- 
ferential equations used to model superfluids and BEC. The description of 
the statistical properties of the solution for large times when a quasi-static 
regime is reached has been developped for a finite number of modes trunca- 
tion. Seeking a stationary probability density function we were led to consider 
the maximum entropy density where all the mass is contained in the coherent 
structure minimizing the Hamiltonian in the continuum limit. Thus, the large 
scale solutions of the dynamics for a finite number of modes are found by this 
approach to behave like this coherent structure immerged in a sea of energy 
equipartitioned fluctuations. Long time numerical simulations show then very 
good agreement with the predictions both for a weakly focusing NLS equation 
and for a ID version of the BEC model. The road towards these statistically 
stationary states gives a consistent scenario of the continuum limit of the 
dynamics. It corresponds to a quasistatic dynamics where the number of fluc- 
tuation modes at equilibrium increase with time. The kinetics of the advance 
towards smaller and smaller scales has been studied for a particular case of 
NLS dynamics and remains to be fully explored in general. 
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